## Load libraries
library(survival)
library(plyr)
library(car)
library(memisc)
library(zoo)
library(gdata)

options(signif.symbols=c("***"=.001,"**"=.01,"*"=.05,"\\dagger"=0.1))

dat <- read.csv("data_used_in_modelling_20130819.csv",header=T)

## Do some factoring
dat$Type<-factor(dat$Type)
dat$Type<-relevel(dat$Type,"Advisory NDPB")
dat$Type2<-car::recode(dat$Type,
	"'Non-ministerial department'='Executive NDPB'")

dat$Function <- relevel(dat$Function,"Representative")

dat$Policy.char<-factor(dat$Policy.char)
dat$Policy.char<-relevel(dat$Policy.char,"Agriculture")


## Full model
mod.cox.full.lag1<-coxph(Surv(start, stop, 1-alive)~
    Function+Policy.char+
	Economic.+I(Function=="Regulation"):Economic.+
	BudgetDelta.l1+#ChangeGov.l1+
	BudgetPos.l1+
	BudgetPos.l1:Debt.l1+Debt.l1+	
	Type+
	exec:start+pubcomp:start+Economic.:start+Debt.l1:start+Debt.l1:I(start^2),
	data=dat)

mod.cox.reduced.lag1 <- coxph(Surv(start, stop, 1-alive)~
    Function+#Policy.char+
	Economic.+I(Function=="Regulation"):Economic.+
	BudgetDelta.l1+#ChangeGov.l1+
	BudgetPos.l1+
	BudgetPos.l1:Debt.l1+Debt.l1+	
	Type+
	exec:start+pubcomp:start+Economic.:start+Debt.l1:start+Debt.l1:I(start^2),
	data=dat)

source("mtable-ext.R")
my.out<-mtable(`Baseline model`=mod.cox.full.lag1,`Reduced model`=mod.cox.reduced.lag1)
my.out<-relabel(my.out,
	"ChangeParty"="Change of Party",
	"ChangeParty.l1"="Change of Party, one-year lag",
	"ChangeParty.l2"="Change of Party, two-year lag",
	"ChangeGov"="Change of government",
	"ChangeGov.l1"="Cumulative changes of government, one-year lag",
	"ChangeGov.l2"="Change of government, two-year lag",
	"CMPDelta"="Change of govt. pos'n., L-R scale",
	"as.factor(COFOG.1): 2"="Defence",
	"as.factor(COFOG.1): 3"="Public order and safety",
	"as.factor(COFOG.1): 4"="Economic affairs",
	"as.factor(COFOG.1): 5"="Environmental protection",
	"as.factor(COFOG.1): 6"="Housing and community",
	"as.factor(COFOG.1): 7"="Health",
	"as.factor(COFOG.1): 8"="Recreation, culture, religion",
	"as.factor(COFOG.1): 9"="Education",
	"as.factor(COFOG.1): 10"="Social protection",
	"Elec.Year"="Election year",
	"Elec.Year.l1"="Election year, one-year lag",
	"CMPDelta.l1"="Abs. change in govt. pos'n., L-R scale, one-year lag",
	"CMPDelta.l2"="Abs. change in govt. pos'n., L-R scale, two-year lag",
	"CMPPos"="govt. pos'n., L-R scale",
	"BudgetPos"="govt. pos'n., L-R scale",
	"BudgetDelta"="Abs. change in govt. pos'n., L-R scale",
	"BudgetDelta.l1"="Abs. change in govt. pos'n., L-R scale, one-year lag",
	"BudgetDelta.l2"="Abs. change in govt. pos'n., L-R scale, two-year lag",
	"BudgetPos.l1"="Govt. pos'n., L-R scale, one-year lag",
	"BudgetPos.l2"="Govt. pos'n., L-R scale, two-year lag",
	"CMPPos.l1"="govt. pos'n., L-R scale, one-year lag",
	"CMPPos.l2"="govt. pos'n., L-R scale, two-year lag",
	"Economic."="NDPB has economic role",
	"Debt"="Debt as percentage of GDP",
	"Debt.l1"="Debt as percentage of GDP, one-year lag",
	"Debt.l2"="Debt as percentage of GDP, two-year lag",
	"BudgetPos x Debt"="Govt L-R pos. * Debt",
	"exec x start"="Executive NDPB * time",
	"start x pubcomp"="Public company * time",
	"Economic. x start"="NDPB has economic role * time",
	"Economic. x I(Function == \"Regulation\")"="NDPB has economic role * Regulation",
	"Debt.l1 x start"="lagged Debt * time",
	"Debt.l1 x I(start^2)"="lagged Debt * time^2",
	"BudgetPos.l1 x Debt.l1"="Govt L-R pos. * Debt, one-year lag",
	"BudgetPos.l2 x Debt.l2"="Govt L-R pos. * Debt, two-year lag",
	summary.stats=c("AIC", "N","Rsq","BIC")
)

capture.output(toLatex(my.out),file="diagnostics/my_models2b.tex")
system("sed -i -e '/Policy.char/,+3d' diagnostics/my_models2b.tex ")

## Now do graphs for interaction terms
my.yhats<-predict(mod.cox.reduced.lag1,type="expected",reference="sample")
my.yhats2<-predict(mod.cox.reduced.lag1,type="lp",reference="sample")
my.yhats3<-predict(mod.cox.reduced.lag1,type="risk",reference="sample")

new.df<- data.frame(Type="Advisory NDPB",exec=0,
	start=mean(dat$start,na.rm=T),
	Function="Representative",Economic.=0,
	ChangeGov.l1=mean(dat$ChangeGov.l1,na.rm=T),
	BudgetDelta.l1=mean(dat$BudgetDelta.l1,na.rm=T),
	BudgetPos.l1=mean(dat$BudgetPos.l1,na.rm=T),Debt.l1=quantile(dat$Debt.l1),
	pubcomp=0)
## plot(predict(mod.cox.reduced.lag1,new.df,type="risk",reference="sample"))

## baseline over time
baseline.df<- data.frame(start=quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),Type=rep("Advisory NDPB",21),exec=0,
	Function="Representative",Economic.=0,
	BudgetDelta.l1=mean(dat$BudgetDelta.l1,na.rm=T),
	ChangeGov.l1=mean(dat$ChangeGov.l1,na.rm=T),
	BudgetPos.l1=mean(dat$BudgetPos.l1,na.rm=T),Debt.l1=0,
	pubcomp=0)
## exec over time
exec.over.time.df<- data.frame(start=quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),Type=rep("Executive NDPB",21),exec=1,
	Function="Representative",Economic.=0,
	BudgetDelta.l1=mean(dat$BudgetDelta.l1,na.rm=T),
	ChangeGov.l1=mean(dat$ChangeGov.l1,na.rm=T),
	BudgetPos.l1=mean(dat$BudgetPos.l1,na.rm=T),Debt.l1=0,
	pubcomp=0)
## pubcomp over time
pubcomp.over.time.df<- data.frame(start=quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),Type=rep("Public companies",21),exec=0,
	Function="Representative",Economic.=0,
	BudgetDelta.l1=mean(dat$BudgetDelta.l1,na.rm=T),
	ChangeGov.l1=mean(dat$ChangeGov.l1,na.rm=T),
	BudgetPos.l1=mean(dat$BudgetPos.l1,na.rm=T),Debt.l1=0,
	pubcomp=1)

baseline.preds<-predict(mod.cox.reduced.lag1,baseline.df,type="risk",reference="sample")
exec.over.time.preds<-predict(mod.cox.reduced.lag1,exec.over.time.df,type="risk",reference="sample")
pubcomp.over.time.preds<-predict(mod.cox.reduced.lag1,pubcomp.over.time.df,type="risk",reference="sample")

pdf(file="diagnostics/interaction_effects.pdf")
plot(quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),baseline.preds,type="b",pch=17,lty=1,xlab="Time",ylab="Risk ratio",ylim=c(0,6),xlim=c(0,25))
lines(quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),exec.over.time.preds,pch=15,lty=2,type="b")
lines(quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),pubcomp.over.time.preds,pch=19,lty=3,type="b")
legend("topright",lty=1:3,pch=c(17,15,19),legend=c("Baseline","Exec","Pub.comp"),bty="n")
dev.off()

## Econ role
econ.over.time.df1<- data.frame(start=quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),Type=rep("Advisory NDPB",21),exec=0,
	Function="Representative",Economic.=0,
	BudgetDelta.l1=mean(dat$BudgetDelta.l1,na.rm=T),
	ChangeGov.l1=mean(dat$ChangeGov.l1,na.rm=T),
	BudgetPos.l1=mean(dat$BudgetPos.l1,na.rm=T),Debt.l1=0,
	pubcomp=0)
econ.over.time.df1$Function<-factor(econ.over.time.df1$Function,levels=c("Representative","Regulation"))
econ.over.time.df2<- data.frame(start=quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),Type=rep("Advisory NDPB",21),exec=0,
	Function="Representative",Economic.=1,
	BudgetDelta.l1=mean(dat$BudgetDelta.l1,na.rm=T),
	ChangeGov.l1=mean(dat$ChangeGov.l1,na.rm=T),
	BudgetPos.l1=mean(dat$BudgetPos.l1,na.rm=T),Debt.l1=0,
	pubcomp=0)
econ.over.time.df2$Function<-factor(econ.over.time.df2$Function,levels=c("Representative","Regulation"))
econ.over.time.df3<- data.frame(start=quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),Type=rep("Advisory NDPB",21),exec=0,
	Function="Regulation",Economic.=1,
	BudgetDelta.l1=mean(dat$BudgetDelta.l1,na.rm=T),
	ChangeGov.l1=mean(dat$ChangeGov.l1,na.rm=T),
	BudgetPos.l1=mean(dat$BudgetPos.l1,na.rm=T),Debt.l1=0,
	pubcomp=0)
econ.over.time.df3$Function<-factor(econ.over.time.df3$Function,levels=c("Representative","Regulation"))

econ.over.time.preds1<-predict(mod.cox.reduced.lag1,econ.over.time.df1,type="risk",reference="sample")
econ.over.time.preds2<-predict(mod.cox.reduced.lag1,econ.over.time.df2,type="risk",reference="sample")
econ.over.time.preds3<-predict(mod.cox.reduced.lag1,econ.over.time.df3,type="risk",reference="sample")


pdf(file="diagnostics/econ_effects.pdf")
plot(quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),econ.over.time.preds1,type="b",pch=17,lty=1,xlab="Time",ylab="Risk ratio",ylim=c(0,6),xlim=c(0,25))
lines(quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),econ.over.time.preds2,pch=15,lty=2,type="b")
lines(quantile(dat$start,na.rm=T,probs=seq(0,1,.05)),econ.over.time.preds3,pch=19,lty=3,type="b")
legend("topleft",lty=1:3,pch=c(17,15,19),legend=c("Baseline","Representative Economic","Regulatory Economic"),bty="n")
dev.off()


## Second lot
## partisanship over debt
baseline.df<- data.frame(start=mean(dat$start,na.rm=T),Type="Advisory NDPB",exec=0,
	Function="Representative",Economic.=0,
	BudgetDelta.l1=mean(dat$BudgetDelta.l1,na.rm=T),
	ChangeGov.l1=mean(dat$ChangeGov.l1,na.rm=T),
	BudgetPos.l1=mean(dat$BudgetPos.l1,na.rm=T),Debt.l1=quantile(dat$Debt.l1,na.rm=T,probs=seq(0,1,.05)),
	pubcomp=0)

## debt over time, left-wing government
left.over.debt.df<- data.frame(start=mean(dat$start,na.rm=T),Type=rep("Advisory NDPB",21),exec=0,
	Function="Representative",Economic.=0,
	BudgetDelta.l1=mean(dat$BudgetDelta.l1,na.rm=T),
	ChangeGov.l1=mean(dat$ChangeGov.l1,na.rm=T),
	BudgetPos.l1=mean(dat$BudgetPos.l1,na.rm=T)-sd(dat$BudgetPos.l1,na.rm=T),Debt.l1=quantile(dat$Debt.l1,na.rm=T,probs=seq(0,1,.05)),
	pubcomp=0)

## debt over time, right-wing government
right.over.debt.df<- data.frame(start=mean(dat$start,na.rm=T),Type=rep("Advisory NDPB",21),exec=0,
	Function="Representative",Economic.=0,
	BudgetDelta.l1=mean(dat$BudgetDelta.l1,na.rm=T),
	ChangeGov.l1=mean(dat$ChangeGov.l1,na.rm=T),
	BudgetPos.l1=mean(dat$BudgetPos.l1,na.rm=T)+sd(dat$BudgetPos.l1,na.rm=T),Debt.l1=quantile(dat$Debt.l1,na.rm=T,probs=seq(0,1,.05)),
	pubcomp=0)


baseline.preds<-predict(mod.cox.reduced.lag1,baseline.df,type="risk",reference="sample")
left.preds<-predict(mod.cox.reduced.lag1,left.over.debt.df,type="risk",reference="sample")
right.preds<-predict(mod.cox.reduced.lag1,right.over.debt.df,type="risk",reference="sample")

pdf("diagnostics/debt_over_partisanship.pdf")
plot(quantile(dat$Debt.l1,na.rm=T,probs=seq(0,1,.05)),baseline.preds,
	type="b",pch=17,lty=1,
	xlab=expression(paste(plain(Debt), " ", (sigma), sep="")),
	ylab="Risk ratio",ylim=c(1,2.5))
lines(quantile(dat$Debt.l1,na.rm=T,probs=seq(0,1,.05)),left.preds,pch=15,lty=2,type="b")
lines(quantile(dat$Debt.l1,na.rm=T,probs=seq(0,1,.05)),right.preds,pch=19,lty=3,type="b")
legend("topleft",lty=1:3,pch=c(17,15,19),
	legend=c("Baseline",expression(paste(plain(Left)," (",-1,sigma,")",sep=" ")),
		expression(paste(plain(Right)," (",+1,sigma,")",sep=" "))),
		bty="n")
dev.off()



## Do some plotting

lp<-by(mod.cox.reduced.lag1$linear.predictors,dat$Name[-c(mod.cox.reduced.lag1$na.action)],mean,na.rm=T)
lp<-data.frame(Name=names(lp),lp=as.vector(lp))
dat<-merge(dat,lp,all.x=T,all.y=F)

topslice.df<-ddply(dat,.(Name),transform,meanlp=mean(lp))
topslice.df<-ddply(topslice.df,.(Name),function(df){df[1,]})

est<-data.frame(Year=names(table(topslice.df$YearEstablished)),Established=as.vector(table(topslice.df$YearEstablished)))
killed<-data.frame(Year=names(table(topslice.df$YearKilled)),Terminated=as.vector(table(topslice.df$YearKilled)))
plot.df<-merge(est,killed)

plot.mat<-t(as.matrix(plot.df[,2:3]))
colnames(plot.mat)<-plot.df$Year

pdf(file="diagnostics/barplot.pdf",width=7,height=5)
barplot(plot.mat,beside=TRUE,las=2,
	legend.text = c("Established", "Terminated"),
             args.legend = list(x = "topright",bty="n",cex=1.2))
dev.off()

